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Plasmon-enhanced Raman scattering can push single-molecule vibrational spectroscopy beyond 
a regime addressable by classical electrodynamics. We employ a quantum electrodynamics (QED) 
description of the coherent interaction of plasmons and molecular vibrations that reveal the emer¬ 
gence of nonlinearities in the inelastic response of the system. For realistic situations, we predict the 
onset of phonon-stimulated Raman scattering and a counterintuitive dependence of the anti-Stokes 
emission on the frequency of excitation. We further show that this QED framework opens a venue 
to analyze the correlations of photons emitted from a plasmonic cavity. 


Surface Enhanced Raman Scattering (SERS) is a spec¬ 
troscopic technique in which the inelastic scattering from 
a molecule is increased by placing it in a hotspot of 
a plasmonic cavity, where the electric fields associated 
with the incident and the scattered photons are strongly 
enhanced (see the schematic in Fig. 1(a)).® The dif¬ 
ference between the energy of those two photons pro¬ 
vides a fingerprint of the molecule, i.e., detailed chem¬ 
ical information about its vibrational structure. Since 
the initia l ob servation of Raman scattering from single 
molecules,®^ the use of a variety of plasmonic struc¬ 
tures that act as effective optical nanoantennas over 
the last decades has allowed a tremendous advance of 
this molecular spectroscop^. Metallic particles such 
as nanoshell^^, nanoringsP, nanorod^, nanowire^, or 
nano star^, as well as plasmonic nanogap structures 
formed in particle dimerJS®^, nanoparticle-on-a-mirror 
morphologie^l^i^, or nanoclusterJi^, are among the vari¬ 
ety of structures that offer huge and controllable enhance¬ 
ments of the field intensity in their hotspots, boo sting the 
inherently weak Raman scattering intensity,and ulti¬ 
mately enabling the chemical identification and imaging 
of particular vibrational modes of a molecule with sub¬ 
nanometer resolution.!^ These results suggest that some 
experiments might have reached the regime where the 
quantum-mechanical nature of both the molecular vibra¬ 
tions and the plasmonic cavity emerges,!^ and call for 
an adequate theoretical description that goes beyond the 
classical treat ment of the electric fields produced in plas¬ 
monic cavities.®22ll3l 

In this work we address the underlying quantum- 
mechanical nature of Raman scattering processes by 
quantizing as bosonic excitations both the vibrations of 
the molecule and the electromagnetic field of a plasmonic 
cavity. The description of the vibrations through bosonic 
operators can be justified by considering the harmonic 
approximation to the energy landscape of the molecule 
along a generalized atomic coordinate (Fig. 1(b)) , such 
as the length of a molecular bond, e.g., G=0.!2il^ These 
vibrations interact with the cavity photons through a 
nonlinear Hamiltonian, reminiscent of that found in op¬ 
tomechanical systems.!^ In this description, the large 


enhancement of the Raman scattering from a molecule 
in the plasmonic cavity occurs thanks to the significant 
shrinking of the effective mode volume of a single photon. 
On the other hand, such plasmonic systems experience 
strong Ohmic losses - that is, they are characterized by 
very low quality factors (typically Q < 10) - placing these 
systems deep in the bad cavity regime (see the sketch of 
a plasmon mode in Fig. 1(c)). This system was recently 
addressed using a linearization scheme based on the clas¬ 
sical Langevin eq uation for the displacements of plasmon 
and phonon fields 

In this work we complement this approach by employ¬ 
ing an exact numerical solution of the quantum mechan¬ 
ical dynamics of the nonresonant Raman scattering pro¬ 
cess in a plasmonic cavity that fully describes the buildup 
of incoherent population of vibrations in the molecule, 
addressing the effect of phonon stimulated Stokes scat¬ 
tering, and furthermore, allows to account for the quan¬ 
tum and classical covariances between the two bosonic 
fields. This formalism thus predicts a variety of nonlinear 
quantum effects in the scattered signal that might have 
already been revealed in different Raman measurements, 
and opens the door for specific design of experimental 
configurations that can test, and eventually control, the 
underlying coherences within SERS. 


I. RESULTS AND DISCUSSION 

A. Theoretical Framework 

The interaction between plasmons and vibrations that 
governs the dynamics of a SERS process can be prop¬ 
erly addressed by establishing the Hamiltonian of the 
system. To derive this Hamiltonian, let us consider the 
vibrational states of the ground electronic state of the 
molecule. In a simplified one-dimensional model the po¬ 
tential landscape of this state can be approximated as 
a displaced harmonic potential. Therefore, the vibra¬ 
tional levels separated by phonon frequency can be 
quantized using the creation and annihilation operators 
b and U. We can thus define a linear polarizability of the 
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FIG. 1 . (a) Schematic of a setup for observing the Raman 
scattering from a molecule placed in a plasmonic cavity, (b) 
Schematic of the two-photon nonresonant Stokes scattering 
between two vibrational states of the molecule (n = 0 —>■ 1) 
mediated by a virtual state |u). A harmonic potential (solid 
lines) approximates the energy landscape of the ground elec¬ 
tronic level (dashed lines), (c) Typical plasmonic cavity mode 
centered at the frequency Wc in the visible range, with the 
width associated with the decay rate k related to lJc through 
the quality factor Q=ujcli^- 


molecule along this coordinate as = R^Q^ib + 
where is the element of the Raman tensor and Qj) is 
the zero-point amplitude of the vibrations. The molecule 
is coupled to the quantized field of the cavity mode 
with resonant frequency lOc and effective volume Veg, 
which can be expressed by the plasmon annihilation (d) 

and creation (d^ operators as E = 'i\J 2 ^° — d^, 
where e is the permittivity of the medium. Thus the in¬ 
duced Raman dipole pa = will be interacting with 
the cavity field E, yielding the interaction Hamiltonian 
Hi = —prE, which, after dropping the fast rotating 
terms oc (d^^ and (d)^, and redefining the equilibrium 
position of the vibrations, yields the non-linear interac¬ 
tion Hamiltonian {H = 1): Hj = —gd^d{w b), with 
the coupling coefficient g = R^Q^oJe/i^Veg). This def¬ 
inition can be related to the resonant emitter-plasmon 
coupling parameter of the Jaynes-Cummings Hamilto¬ 
nian gjc and the Purcell factor Ep oc QjVeg (see the 
discussion in Supporting Information). Interestingly, the 
interaction Hamiltonian Hi is identical to the one used 
in optomechanical systems, in which the quantized oscil¬ 
lations of a c avity mirror modify the resonance frequency 
of the cavity.l22E5l 

It should be noted that this approach is limited to 
the off-resonant Raman scattering, for which the vir¬ 
tual state mediating the Raman transition (Fig. 1(b)) is 
strongly detuned from an excited electronic state. Also, 
our Hamiltonian does not consider Raman processes that 
involve the effect of the simultaneous excitation or emis¬ 
sion of two or more phonons, which may become impor¬ 


tant for very intense lasers. 

Furthermore, our model assumes that the plasmonic 
system can be approximated as a single cavity mode, and 
described through the canonical quantization scheme. 
Although this approach should yield an accurate descrip¬ 
tion for the many of plasmonic SERS setups, a more exact 
description could be offered by a rigorous quantization 
method which tak es int o account the exact Green’s func¬ 
tion of the systerrpS1291 and would allow consideration of 
complex plasmonic responses involving many modes. 

The coherent evolution of the system, including an ad¬ 
ditional driving term (laser) at frequency w/ is then given 
by the full Hamiltonian 

H = + —50^0(6’^ -1-6) -|-iH(a’^e““'‘ — de*“'*). 

( 1 ) 

Throughout the paper will be referred to as a pumping 
power proportional to the power density of the input laser 
and the intrinsic parameters of the cavity mode. 

The master equation for the dynamics of the density 
matrix of the system p in the absence of pure dephasing^Sl 
reads 

_ T Ai , r 1 , + l)7m„ r , , 

dtp — t[p, H]-{--T>a[p]-\ - - --^ [p]i 

( 2 ) 

where the last three terms on the righ t-han d side are 
Lindblad-Kossakowski terms defined aPi^^ ^o[p] = 
20pO^ — O^Op — pO^O, describing the decay of the 
plasmons (27a), phonons (27j,), and the thermal pumping 
(27g^) of the phonons by the environment at temperature 
T, with the thermal population = (^ff^-rn/ksT _ i)-i. 

In the following, we take the phonon energy tOm = 0.1 
eV and the phonon decay rate 7m = 1 meV,l^ as rep¬ 
resentative values that characterize vibrational finger¬ 
prints, and consider a typical plasmonic cavity with qual¬ 
ity factor of Q = 10 resonant at oje = 2.5 eV, decaying 
with rate k = 0.25 eV (see Supp. Info.). The coupling 
parameter is taken as g = 1 meV in accordance with 
the reported characteristics of state-of-the-art resonant 
plasmon-emitter systems (see Supp. Info.). 

To analyze the dynamics of the system, we follow 
three complementary approaches. First, we solve di¬ 
rectly the master equation in Eq. ([^ by means of nu¬ 
merical calculations,!^ (taking advantage of low excita¬ 
tion numbers to truncate the Hilbert space). We comple¬ 
ment this numerical solution by two other approximate 
treatments that allow to obtain analytical results, use¬ 
ful to identify the influence of the different parameters 
in the Raman signal and thus provide simpler expres¬ 
sions to interpret experimental results. In the first an¬ 
alytical approximation, we linearize the Hamiltonian to 
a purely quadratic form which applies to the weak cou¬ 
pling {g k) regime we are interested in, and allows 
for the analytical treatment (see Supporting Information 
for details). We note that when solving the linearized 
Hamiltonian, we do not map the quantum Langevin 
equation to the classical dynamics equations, as is often 
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done in the analysis of optomechanical systems!^ Con¬ 
sequently, this approach provides a complete character¬ 
ization of the classical and quantum correlations within 
the system. Finally, to develop the second analytical 
approach, we apply the quantum noise approach^ in 
which the vibrations of the molecule are coupled to a 
bath of fluctuating population of the cavity plasmons.^^ 
This last formalism allows us to describe the effects of 
the static and dynamical backaction, which arise when 
the vibrations of the molecule modify the resonant fre¬ 
quency of the cavity and, the refore , the population term 
in the Hamiltonian (Eq. The static backaction 

describes the constant, coherent displacement of the state 
of the vibrations, and can be largely neglected in the 
regime of parameters of interest (see the Methods sec¬ 
tion for a further discussion). On the other hand, the 
dynamical backaction (DBA), observed when the inci¬ 
dent laser is detuned from the cavity, leads to the heat¬ 
ing or cooling of the vibrations.^ The quantum noise 
approach captures those effects, at the price of losing the 
inform ation about coherences in the molecule-plasmon 
system .1231311 The spectra of emission from the cavity are 
calculated following Glauber’s photodetection theory,!^ 
as /(w) = Qfdet<S'(a;), where the frequency-independent 
parameter Odet describes the properties of the detection 
systenJ22l and S{uj) = dte-“‘(at(t)a(0)),, in the 

steady state (ss). The two-time correlator is calculated 
by applying the quantum regression theorem (QRT).l^ 

Let us first focus on the numerical solution. Two spec¬ 
tra of emission from the cavity, calculated for the weak 
(H^ = 10“^ eV^) and strong illuminations (H^ = 0.5 eV^) 
and for the laser tuned to the cavity (A = Wc — w; = 0), 
are shown in Fig. 2(a) and (b), respectively. In the in¬ 
set of Fig. 2(a) we zoom in on the anti-Stokes emission 
calculated for the environment at T = 0 K (dashed line) 
and T = 300 K (solid line; nf'iT = 300 K) « 0.02). 
The difference between these plots illustrates the ef¬ 
fect of ther mal pu mping of the vibrational levels by the 
environment .ESESl 

To further explore how the signal evolves with the 
pumping power, we plot in Fig. 2(c) the maxima of 
the Stokes (blue solid line; in the region of interest this 
function is independent of the temperature T) and anti- 
Stokes (orange dashed line for T = 0 K and orange solid 
line for T — 300 K) emission lines for the increasing 
In the weak pumping regime (H^ < 10“^ eV^) for 
nonzero temperature (solid lines) both the Stokes S{ujs) 
and anti-Stokes S{ujas) emission depend linearly on 
indicating that the anti-Stokes transition originates from 
the thermally excited vibrational state. For higher driv¬ 
ing powers (10“^ eV^ ^ ^ 0-5 eV^) the anti-Stokes 

intensities become independent of the temperature, as 
the phonons are provided primarily by the Stokes tran¬ 
sitions {vibrational pumping)W^ The transition between 
the thermal and the vibrational pumping of phonon^^ 
is illustrated in Fig. 2(d), where we plot the populations 
of the phonons (green line) and plasmons (red line) for 
T = 0 K (dashes lines) and T = 300 K (solid lines). 
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FIG. 2. Dependence of the Raman scattering on the exci¬ 
tation power and temperature: (a,b) emission spectra S{u)) 
of the molecule in the plasmonic cavity for pumping (a) 
= 10“^ eV'^ and (b) = 0.5 eV^ at T = 0 K (dashed 

lines) and T = 300 K (solid lines); (c) emission intensities of 
the Stokes {S{uJs), blue lines, calculated at T = 300 K) and 
anti-Stokes photons {S{uJas), orange dashed and solid lines 
for T = 0 K and T = 300 K, respectively) as a function of the 
pumping power (d) populations of plasmons (red line) 
and phonons (green lines) in the steady state for T = 0 K 
(dashed lines) or T = 300 K (solid lines). Vertical dashed 
lines in (c,d) indicate the pumping powers for plots in (a) and 
(b). All the cases assume A = 0. 


Finally, for the highest considered pumping powers 
(H^ > 0.5 eV^) the Stokes intensity S{ojs) visibly sur¬ 
passes the expected linear dependence on (marked 
with dotted gray line in the top-right corner of Fig. 2(c)). 
To address this effect, we develop the first of the afore¬ 
mentioned analytical approaches, by considering the lin¬ 
earized Hamiltonian and employing the QRT which, for 
the particular case of A = 0,7m “C uJm,K and environ¬ 
ment temperature T, yields the following expressions for 
the Stokes and anti-Stokes photons: 

5(a;s) = + + —V (3) 

7m \ 7m / 

and 

S{ 0 Jas) = ^S2H2 (nt + s^n^—) , (4) 

'Tm \ 'Im J 

with S 2 ~ [4:g/{K\K — 2iuj^\)]^ (see Supp. Info, for the 
complete derivation). The first term in the brackets of 
Eq. ([^ can be recognized as the conventional two-photon 
cavity-assisted Stokes transition, linearly dependent on 
As we show in the Supp. Info., by considering the en¬ 
hancement K = \E/Eq\ of the electric component of the 
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incident illumination Eq produced by the plasmonic cav¬ 
ity, and relating it through the re ciprocity theorem to the 
Purcell factor of the cavity, 1^2111] this term in Eq. (|^ al¬ 
lows us to retrieve the expected dependence of the Stokes 
emission S{ujs) oc K'^. 

The sum of the second and the third terms in brack¬ 
ets in Eq. ([|), or the two terms in brackets in Eq. Q, 
represents the incoherent population of the phonon mode 
= (S'l'S)ss — {b^)ss{b)ss) arising from (i) the ther¬ 
mal pumping and (ii) the coupling to the plasmon cav¬ 
ity, respectively. These two terms together describe a 
process of phonon-stimulated Raman scattering, in which 
the population of vibrations enhances the rate of Stokes 
scattering. Therefore, the last term in Eq. ([^ explains 
the nonlinearity of the Stokes signal for large We 
note that phonon-stimulated Raman scattering has been 
reported in experiments on ensembles of Raman-active 
centers, e.g., hydrogen gas (see Refs. [142143] and refer¬ 
ences therein). However, to our knowledge, the criterion 
for the onset of this effect for single scatterers in a cavity 
(see the Supporting Information for an explicit formula¬ 
tion) has not been reported in the literature. 

To further explore the effects of the thermal and vi¬ 
brational pumping of phonons on Raman scattering, we 
consider the dependence of the Raman scattering on the 
detuning A = Wc —w; of the incident laser. In the typical 
classical models of SER^^ the dependence of the Stokes 
( 5 'ciass(^^; anti-Stokes (S"='^""(a;as; A)) emis¬ 

sion is determined by the enhancement of the electric 
field of both the incoming {\E{uJi)/Eo{uJi)\‘^) and outgo¬ 
ing {\E{oJs/as)/Eo{uJs/as)?) photons at the position of 
the molecule: 


^^'-^(a;s/a5;A)cx4/,s 


E{u;i) E{u}s/as) 
Eo{uJi) Eo{u!s/as) 


(5) 


Assuming that the enhancement is given by a Lorentzian 
profile centered on the cavity resonance ujc with width k 
(gray curves in the bottom panels of Fig. 3), we expect 
that the calculated Raman emission spectra S{u!s/aSj A) 
will depend on the laser frequency as depicted with the 
blue dashed curves in the bottom panels of Fig. 3. In 
particular, the Stokes signal should be strongest for the 
incident laser blue-tuned from the cavity. This general 
result for the Stokes scattering is confirmed by our nu¬ 
merical calculations performed for various pumping pow¬ 
ers (H^ = 10“^ eV^ to 0.25 eV^, from red to green lines) 
and different temperatures (T = 0 K as dashed lines and 
300 K as solid lines), as shown in the upper panels of 
Fig. 3(a). 

However, there are also marked differences: we note 
that when the laser is blue-detuned from the cavity res¬ 
onances (A < 0), the profile S{u}s; A) becomes narrower 
as we increase the pumping intensity. This behavior re¬ 
sults from the buildup of the incoherent population of 
phonons due to the vibrational pumping and to the dy¬ 
namical backaction.l^ The latter phenomenon, triggered 
only when the laser is blue-detuned from the cavity res¬ 
onance (A < 0), has been observed in optomechanical 


systems,!^ and was recently argued to occur in Raman 
scattering from molecules.!^ We discuss this concept in 
more detail in the following paragraphs. 

Interestingly, the classical Eq. fails to explain the 
dependence of the anti-Stokes scattering obtained in the 
full quantum calculations (Fig. 3(b)). For the weakest 
driving powers (red lines, = 10“^ eV^), the S{u}as] A) 
intensity is found to be the largest for the laser frequency 
uji on resonance with the cavity lUc in the absence of ther¬ 
mal pumping of phonons (T = 0); however this max¬ 
imum appears as red-tuned from the cavity resonance 
when the thermal pumping of phonons is produced (at 
T = 300 K). For stronger driving powers, the intensity 
plots for T = 0 K and 300 K start to merge and peak 
at increasingly blue-shifted frequencies, notably crossing 
the cavity resonance. As for the supra-linear increase 
of the Stokes scattering, this surprising property stems 
from the laser pumping of the vibrational levels. The 
classical Eq. ([^ for the anti-Stokes A) inten¬ 

sity does not account for the origin of the vibrations of 
the molecule, and therefore can only be applied when 
these are provided by the heated reservoir (it should be 
noted, however, that a suitable correction to ^(was; A) 
introducing the vibrational pumping has been proposed 
by Kneipp et 

To understand the dependencies of the Stokes and anti- 
Stokes emission on the detuning presented in Fig. 3(a) 
and (b), we need to consider the vibrational state of the 
molecule in detail. Since the linearized Hamiltonian can¬ 
not be easily solved analytically for an arbitrary detuning 
A, we develop an alternative analysis of the dynam ics of 
the system based on the quantum noise approach,l23EIl 
which also allows to obtain analytical expressions for the 
Stokes and anti-Stokes photon emission. According to 
this approach, and similarly to the weak coupling model 
of a two-level system in a cavity, the rates of relaxation 
and excitation of vibrations are considered to be pro¬ 
portional to the noise spectrum determined by the char¬ 
acteristics of the plasmonic cavity. This coupling effec¬ 
tively modifies the decay rate of the molecular vibrations 
to 7m + 7opt) where the optomechanical damping 7opt is 
dependent on the detuning A and the intensity of the 
incident laser (see Eq. (15)). Furthermore, 7opt describes 
the exchange of energy between the vibrations and fluc¬ 
tuations of the cavity population about its mean value. 
A description of the quantum noise approach with the 
explicit expression of 7 opt [Eq. ([T§] is given in the Meth¬ 
ods section. As a result, the intensities of the Stokes 
and anti-Stokes emissions for general detuning are given 
within the quantum noise approach by: 


S{ujs) 




^incoh 


+ 1 


(A -|- LUjn)^ + (^/2)^ 7opt + Ira 


and 


S{uJas) OC 




Tincoh 


( 6 ) 


(7) 


(A - UJrnY + (k/2)^ 7opt + 7m ’ 
respectively, where denotes the incoherent phonon 
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FIG. 3. Dependence of the (a) Stokes and (b) anti-Stokes emission and (c) phonon population on the frequency of the incident 
laser wi. In the two left-most top panels we show the numerically calculated intensities at the (a) Stokes S{ujs', A = uJc — oji) 
and (b) anti-Stokes S{LOaS', A) frequencies for the pumping power from = 10“^ eV (red lines) to = 0.25 eV (green lines) 
and environment temperature T = 0 K (dashed lines) and T = 300 K (solid lines). In bottom panels in (a, b) the dashed 
blue lines show the predictions of classical models (Eq. ([^) for (a) A) and (b) 5‘^*“°(aJas; A), whereas the cavity 

resonances are shown in (a-c) in gray. In (c) we plot the incoherent phonon population calculated within the quantum 

noise approach [Eq. (141]. Identical results are obtained with the exact numerical calculation. The temperature and intensity 
for each situation are denoted with the same color code as in (a) and (b). 


population, and n™** is the coherent population of the 
cavity (see the Methods section for the exact definition 
of these parameters). The dependence of as a func¬ 
tion of detuning and laser intensity, shown in Fig. 3(c) 
and obtained from application of Eq. (14), helps to un¬ 
derstand the dependences of the Stokes and anti-Stokes 
emission in Figs. 3(a) and (b). These analytical results 
are indeed identical to those obtained with the exact nu¬ 
merical solution. 


For the laser blue-tuned from the plasmonic cavity 
(A < 0), the optomechanical damping 7opt becomes neg¬ 
ative, resulting in an increase of the incoherent phonon 
population as observed in Fig. 3(c). Consequently, 

the Stokes and anti-Stokes emission reveal similar depen¬ 
dence on A, as shown in Figs. 3(a) and (b). The Stokes 
and anti-Stokes peaks’ widths exhibit also a narrowing 
from 7 m to 7m + 7opt^ which, interestingly, is compen¬ 
sated for by the denominators in the above equations, 
ensuring that the integrated intensity of the inelastic 
scattering is only implicitly dependent on 7opt through 
nif,ncoh. Furthermore, we note that in analogy to other 
optomechanical systems, we can define a cooperativity- 
like parameter C = 7opt/7m of the Raman process. For 
the parameters discussed throughout this work and for 
blue-tuned laser, C reaches the minimum value of —0.15 
for A « —O.Swm (see Fig. S2 in the Supp. Info.). For 


C < — 1 (in the case of laser blue-tuned from the cavity 
A < 0), the system would exhibit phonon lasing.^ 

On the other hand, for red-tuned laser (A > 0) we 
enter the so-called cooling regime in optomechanics, in 
which the vibrations can be suppressed below the ther¬ 
mal population However, because we are not in the 
sideband-resolved limit {k Wm), in a typical SERS con¬ 
figuration this effect is suppressed, and does not 

significantly drop below as illustrated by the solid 
lines in Fig. 3(c). 

Last, we focus on an important piece of information 
regarding the characterization of Raman photon emis¬ 
sion that can be obtained from the time- and frequency- 
resolved correlatioiPS®^ between Stokes and anti-Stokes 
photons emitted from the cavity. In the steady state, 
this magnitude can be theoretically calculated through 
intensity-intensity correlations of the filtered output field, 

(• (t)^J,2 ^^(t + T)H^2_r2(i + T)^t.cJi,ri (^)] 0 

p^i<^^,rj(t))((ij,^ p^i^,,r2)(i + t)) 

( 8 ) 


where T and : denote the ti me a nd normal or¬ 
dering operators, respectively.^^^I^ Furthermore, 
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FIG. 4. (a) Two-photon frequency-resolved Stokes-anti- 

Stokes correlators grpp(a;s,tJas) and (b) physical spectra of 
emission = {AJ^ p(0)Ti.j,r(0)) calculated for the tem¬ 

peratures T = 0 K (dashed lines) and T = 300 K (solid 
lines), coupling parameters chosen as multiples of go = k/250: 
g = go (blue lines), g = 2go (orange lines), g = 4go (green 
lines), and laser tuned to the cavity A = 0. 


detected at frequency Wi, within a frequency window 
at time t. For simplicity we consider Fi = r 2 = F 
and place the Lorentzian filters at Wi /2 = ^s/aS- The 
( 2 ) 

photon correlations gj, p (ws, WaS; r = 0), calculated 
using a recently developed methocP^, are plotted in 
Fig. 4(a) as a function of the driving parameter 
for an environment temperature of T = 0 K (dashed 
lines) and T = 300 K (solid lines). The details of the 
method of calculation of the correlation function are 
provided in Supp. Info. The coupling parameters g are 
considered as multiples of go = At/250, namely g = go 
(blue lines), g = 2go (orange lines), g = 4go (green 
lines), with the filter linewidth F = 0.1 meV. As shown 
in Fig. 4(b), for these parameters the physical spectrum 
of emission S'p^^(w) = p^ (0)A^,ri (0)) is formed by 
three peaks: the elastic Rayleigh scattering and inelastic 
Stokes and anti-Stokes contributions (note that S'p^^ does 
not include the explicit dependence on the factor). 
We clearly observe in Fig. 4(a) that for weak coherent 
pumping and in the absence of thermal pumping (dashed 
lines), the system exhibits strong bunching stati stics, 
which is a signature of strongly correlated emission !^^ ! ^^ ! 
The physical origin of the strong correlation is discussed 
in detail in the Supp. Info. In short, in the absence 
of other sources of excitations, the Raman photons are 
emitted by exchanging a single phonon and are therefore 
strongly correlated. This process should exhibit a 
5 p^p(a; 5 , Was) oc 1/S(ujs) oc dependence, as the 
rate of coincidences in the detec tion is determined by 
the anti-Stokes emission The results in 

Fig. 4(a) show this evolution of the correlation for 
large incident power. As the thermal pumping (solid 
lines in Fig. 4(a)) grows, the anti-Stokes photons are 
increasingly created through the absorption of phonons 
which originate from thermal excitation, leading to the 
quenching of the correlation. Finally, the dependence 
of the correlations on the coupling parameter stems 


from the fact that for small g (and weak anti-Stokes 
emission), the filters detect primarily the elastically 
scattered photons (see Fig. 4(b)), leading to a plateau 
in (;p^p(w 5 ,a;as) for small 11. We note these results are 
consistent with those recently reported by Kasperczyk et 
alfM from measurements of the Stokes and anti-Stokes 
pairs emitted from a thin layer of diamond and by 
Parra-Murillo et al!^ from the theoretical analysis of 
the Raman scattering. 

II. CONCLUSIONS 

In summary, we have implemented a fully quantum- 
mechanical model of the inelastic nonresonant Raman 
scattering from a molecule placed in a lossy plasmonic 
cavity. Our analytical and exact numerical calculations 
point to phenomena which go beyond the standard de¬ 
scription of Raman scattering such as (i) the onset of 
stimulated Stokes emission, and (ii) the counterintuitive 
dependence of the anti-Stokes signal on the detuning 
of the incident laser from the cavity resulting from the 
pumping of the molecular vibrations. Furthermore, we 
demonstrate the strong correlations of Stokes and anti- 
Stokes pairs emitted from the cavity and analyze the ori¬ 
gin of this effect. Finally, we note that our model can 
be readily applied to a range of physical problems, i.e., 
investigating the Raman scattering of noncoherent light 
or the quantum correlations induced in the system. 

III. METHODS 

A. Linearization of the optomechanical 
Hamiltonian 

Let us consider the optomechanical Hamiltonian given 
in Eq. Q. If we decouple the cavity from the molecule 
(putting g = 0) the plasmonic system is driven into a 
steady coherent state with amplitude as = H/(| -|- lA). 
The interaction with these steady-state plasmons domi¬ 
nates the dynamics of the system in the weak-coupling 
regime as can be seen by the transformation to a dis¬ 
placed basis (with |as) — |0)) in the plasmonic system. 
By, correspondingly, replacing a —>■ a + as, we can ex¬ 
press the coherent dynamics of the system through the 
Hamiltonian 

H =lS.(A a + — g\as\’^{b+ tA) + 

— g{as(A + a*d)(6 -I- 6^) — geAaih + 6^). (9) 

As we discuss in the Supporting Information in detail, 
in the regime of parameters used throughout this work, 
characterized by weak coupling g, the last nonlinear term 
can be neglected, allowing us to write the simplified 
Hamiltonian 

H' = Aa^ d+uJmii^ b—g\as\'^ {b+b^) — g{ascA +a*a){b+b^), 

( 10 ) 
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which yields a purely quadratic/Gaussian dynamics and 
linear quantum Langevin equations. We will refer to H' 
as a linearized Hamiltonian. 

In this new Hamiltonian, the coherent driving of the 
cavity gives rise {via the nonlinear interaction) to a linear 
couplin g be tween phonons and plasmons (the last term 
in Eq. (lOl) and an effective coherent driving (?|q!sP of 
the phonon mode. To remove this driving term, we dis¬ 
place the phononic operators in the same way as for the 
plasmon, and obtain a new coherent driving of the cav¬ 
ity (and a renormalized detuning), and so successively. 
We can capture all orders of this feedback by defining 
displacements a/ and /3( through the condition that in 
the displaced basis the Hamiltonian does not contain any 
linear {driving) terms. Such displaced Hamiltonian can 
be written as 

H" = A'a^a -I- uimb^b — g\a'^a^ + (a(,)*d](6 -|- S^), (11) 
where displacements a), and /3( are defined by 
, Vl 


k/ 2 + iA' ’ 


( 12 ) 


P's = 




and the renormalized A' = A — 2gYie{fi{). 


(13) 


B. Quantum noise approach 


In the quantum noise approach, the molecular vibra¬ 
tions are coupled to the plasmonic cavity described by the 
quantum noise density of the plasmon population fluctu¬ 
ations a^a defined as Spp{uj) = dte“*(E(t)F(0)), 

where F = cL^d — {a^d) denotes the fluctuations of the 
plasmon number around the equilibrium, and the average 
is taken over the state of the uncoupled cavity {g = 0). 
From Spp we can obtain the effective relaxation and ex¬ 
citation rates of the vibrations. As a result, the molecular 
vibrations are in a displaced thermal state. The displace¬ 
ment given by /3(, and the incoherent population of the 
thermal state is 


incoh 

Ub 


Im th 
“r "Topt 


Topt 
“1“ Topt 


^rad5 


(14) 


The optomechanical damping is defined as 


7opt — g \_Spp{oJm) Spp{ Wm)] — 

1 




(A - Wm)2 -p (k/2)2 


1 


(A -|- Wm)^ + (k/2)2 
(15) 


with n™'' = |a(,p denoting the coherent population of the 
cavity. Here we simplified the F operator by considering 
the displaced cavity operators d —>■ a -I- a), and dropping 
the quadratic term: F « •\/n^(a-|-a'^).I^The first (sec¬ 
ond) term in the brackets of Eq. corresponds to the 
transfer of energy from vibrations to the cavity (from the 
cavity to the vibrations). The more general form of 7opt 
includes the static displacement of vibrations, which can 


be accounted for by replacing the detuning A in Eq. (15) 


by A' defined earlier. However, in the regime of pa¬ 
rameters of interest in a typical Raman configuration, 
this correction is negligible. The second contribution in 
Eq. (14) - the contribution to the population of vibra¬ 


tions originating from the coupling to the optical cavity 
- is proportional to rirad: 


^rad — 


( A -02 + («/ 2)2 


4Awi 


(16) 


Here is the me chani cal frequency corrected for the 
optical spring effect^^^ This correction however can be 
neglected for the parameters in the regime of interest in 
SERS. 

Finally, we note that in the limit of laser tuned to 
the cavity resonance A —)■ 0, there is no DBA mecha¬ 
nism and the optomechanical damping 7opt vanishes. In 
this case, and neglecting the static displacement of vi¬ 
brations (A « A'), we reproduce within the quantum 
noise approach the simple expression for the vibrationally 
pumped incoherent population of vibrations derived ear¬ 
lier from the full quantum-mechanical linearized treat¬ 
ment, thus making both analytical derivations fully con¬ 
sistent: 


lim 

A-i-O 


„th 

Ub 


4ggn0°^ K 
+ 4w^ 7m 


+ — - (17) 

Tm 


Similarly, the expressions for the photon emission inten¬ 
sities S{oJs) and S{u}as) obtained in Eqs. § and 0 with 
application of the quantum noise approach, also match 
those obtained from the previous analytical derivation 
(Eqs. ([^ and @). 


(see the previous section and the Supporting Information 
for details). The second term in Eq. (141 describes the 
effective heating of the vibrations by the incoherent pop¬ 
ulation of the cavity (associated with the fluctuations of 
its population), and dominates the population of phonons 
in the vibrational pumping regime. It should be noted 
that this important contribution to cannot b e de- 

rived by assuming a classical form of the cavity field 
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Appendix A: Derivation of the parameters 
1. Cavity-vibrations coupling 


Here we provide estimates of the coupling parameter g between the Raman transitions in a molecule in vacuum 
and a plasmonic cavity. We show that this parameter can be related to two magnitudes discussed extensively in the 
literature: the Purcell factor Fp and the coupling parameter gjc of the resonant emitter-plasmon Jaynes-Cummings 

Hamiltonian.EMlol 

The coupling parameter g derived in the main body of the paper is given by 

9 = RuQI^, (A1) 

EoPeff 

where R^, is the relevant element of the Raman tensor, Q® is the zero-point amplitude of the vibrations, and the 
fraction on the right-hand side includes the parameters of the cavity mode - its resonance frequency Wc and effective 
volume VeB- Therefore, we can relate it to the Purcell factor of the mode Fp defined as 


Fp 


oJc \ 3£o / 27rc\ ^ 
eoKff/ V Wc / 


arriving at 





(A2) 


(A3) 


Furthermore, we note that in the derivation of the coupling parameter g we have followed a similar path as for 
deriving the resonant coupling gjc between a two-level system and a cavity mode in the Jaynes-Cummings model. For 
the resonant, properly placed and oriented dipolar moment d of the transition of the two-level system, the maximum 
coupling can be expressed as!^ 


\9jc\ = 


huc 
h \ 2eot4ff’ 


(A4) 


where ^/hu}l./{‘2eoVeff) yields the electric field as E = i\J 2^v — a^). Since the Raman processes can be viewed as 

a two-step process mediated by the virtual state, the coupling g is proportional to the square of gjc-, as can be shown 
from Eqs. (A3) and (A4): 


9 = 


R.Ql 




(A5) 


We can therefore find an estimate of the absolute value of the coupling c oeffic ient g by recalling some of the reported 
Purcell factors Fp for cavities resonant in the visible and near-IR regime 


Reference 

Type of cavity 

Fp 

hijjc [eV] 

hK, [eV] 

hg (R6G) [eV] 

Esteban et 

plasmonic patch 

antenna 

5 X 10 ^ 

0.45 

0.025 

2.5 X 10 "®® 

Chen et 

plasmonic particle 

on a dielectric substrate 

8 X 10 ® 

1.5 

0.1 

6 X 10 “® 

Liu et 

plasmon mode 
in nanorod 

5 X 10 ® 

2.2 

0.07 

8 X 10 “® 

Zeng et al¥^ 

plasmonic dimer 

3 X 10 ® 

3.5 

0.2 

6 X 10 “®’ 


TABLE I. Values of the Purcell factors Fp for various types of cavities with resonant energies tiojc and widths Hk. In the 
last column we provide the estimates of the coupling coefficient hg calculated for the coupling with the specific transition of a 
rhodamine 6G molecule (see text below for details). 

The values of the elements of the Raman tensor R^, and the zero-point amplitude Q® vary significantly for different 
molecules and each vibrational degree of freedom. The estimates of these values have been provided from theoretical 
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and experimental studies for numerous molecules, including rhodamine 6 g (R 6 GJP^ or various peptides,^ as well as 
cluster structures of e.g., silicon!^ To provide estimates of the values of g we consider the specihc values of Raman 
activity of rhodamine 6 g. The non-resonant Raman spectra of the R 6 G molecules exhibit vibrational energies in the 
range of one hundred of meV and, in the strong-vibrations end of the parameters spectrum, a Raman activity close 

O 4 

to 5 X 10^ EqA amu“^. In the one-dimensional model of vibrations used throughout this work, the Raman activity is 
equal to R^. Thus, the product where Q® = reaches for t he v ibrational energy hjJm = 0.1 eV the 

value of approximately 3 x 10“^*^ Eq m^. Including this parameter into Eq. (A5), we arrive at the upper estimate of 
the coupling parameter of around hg Ki 6 x 10“^ eV for a molecule placed in a plasmonic dimer nanoantenna.l^ 

We should note that, except for the dark plasmon mode in a nanorod,^ the structures listed above are designed 
to provide large Purcell factors while retaining high quantum efficiency and avoiding quenching of emission from the 
resonant two-level systems. The latter is of limited relevance when designing systems for SERS or TERS (Substrate- 
or Tip-Enhanced Raman Scattering), since the Raman scattering does not suffer from quenching and, therefore, we 
can consider other setups, e.g. metallic dimers with subnanometer gaps and significantly reduced mode volumes. We 
can estimate a lower limit for the mode volume from the recent contribution by Barbry et where the TDDFT 
(Time-Dependent Density Functional Theorem) was employed to analyze the response of a dimer of sodium clusters. 
From the spatial profile of the field localization in the gap, and assuming that all the energy is confined to the gap. 


we get Idn 


10 


-28 


m'^ for the mode energy of hiOc = 3.3 eV. 


To e stim ate the coupling parameter expected in such systems, we can apply the definition of the coupling given 
in Eq. (All. Taking a more conservative value of the mode volume of 14ff = 10“^® m^, and considering the Raman 
activity of rhodamine 6 g, we obtain hg ke \ meV, which we use throughout this manuscript. 


2. Approximations in the derivation of interaction Hamiltonian 

The complete expression for the interaction Hamiltonian Hj = —pnE., obtained by inserting the expressions for the 
Raman dipole pr and the electric field E, is 

Hi = -g jata + ^[-a^ - + 1]| (St + S). (A6) 

Within the rotating wave approximation we remove the first two terms in the square brackets, which yield rapid 
oscillations at frequencies of 2wc 3> uim- Furthermore, we can displace the equilibrium position of the vibrations to 
account for the interaction between vibrations and the vacuum field of the cavity (—|(S -|- S^)). 

Furthermore, we note that by using the free-space Raman tensor R^ which describes the interaction of a molecule 
with an incident planewave in free space, we focus our attention on the electromagnetic enhancement mechanism, 
thus neglecting any contribution from the chemical enhancement. 

Finally, for the sake of simplicity, we neglect the unperturbed component of the linear optical polarizability of the 
molecule. 


3. Coupling of the cavity to the incident light 


The coherent driving parameter of the cavity is defined as* 


len 




(A7) 


where \E^\ is the maximum of the near field scattered by the structure, which can be introduced through the field 
enhancement ability of the plasmonic system as K = |E^|/|Eo|, where |i?o| is the amplitude of the incident coherent 
illumination. The plasmonic cavity can be also characterized by the Purcell factor Fp, the radiative rate enh ancem ent 
r^j/To and the quantum yield 77 of a dipolar emitter placed in such cavity. Through the reciprocity theoremj ^** * ^*^ * and 
for the dipolar mode, these three quantities can be related to K as 


^=K^ = VFP. 

t 0 


(AS) 


Inserting the definition of the Purcell factor (Eq. (A2)) into the second equality, and plugging the resulting expression 
for K into Eq. (A71, we arrive at 


47r 


SgpAg 

2n 


r]K\Eo\. 


(A9) 
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To provide some exemplary values of the parameter we consider an optical plasmonic cavity with the plasmon 
frequency huc = 2.5 eV, quality factor Q = 10 {Hk = 0.25 eV) and a conservative quantum yield rj = 0.01. Thus, for the 
strong laser intensity I = 10^ W/cm^, which yields |i?o| « 8 x 10® V/m, we get the driving parameter Ml « 0.11 eV. 
Nonetheless, stronger values are necessary to observe clearly the effect of the stimulated Raman scattering in the 
absence of the dynamical backaction. Therefore, we should consider either a plasmonic cavity providing a stronger 
enhancement of the incident field, or a stronger illumination source - possibly replacing the continuous illumination 
source with a pulsed laser (capable of enhancing the electric field by an order of magnitude) with the pulse length 
larger than the characteristic time of the buildup of phonon population oc Alternatively, the onset of non¬ 

linearity could be observed with the current continuous illumination systems by considering stronger coupling with 
the molecules, either due to the smaller effective volume of the mode, or the stronger intrinsic Raman activity of the 
probed system {e.g. by replacing the R6G with graphene nanosystem or carbon nanotubes). Interestingly, Roelli et 
cite very similar restrictions for the observation of the parametric instability. 

For the strongest laser intensities used throughout the paper, we approach the regime where and ujc become 
comparable. Nevertheless, even for those strongest pumping coefficients, the pumping term of the Hamiltonian 
(Eq. (2) in the main text) is written in the rotating wave approximation. We expect that the inclusion of the counter¬ 
rotating terms -I- will not change significantly the dynamics of the system, as they will probably 

result mainly in energy shifts that can be corrected with appropriate laser detuning. 


Appendix B: Linearized Hamiltonian and the analytical treatment of the system 

In the following section we present an analytical treatment of the system, based on the quantum regression theorem 
(QRT).I^To readily apply this method, we first consider a linearized version of our Hamiltonian and verify numerically 
that this approximation is well-suited for describing the system in the range of parameters discussed here. We 
then proceed to address the linearized Hamiltonian analytically and present a solution for the case of the incident 
illumination tuned to the plasmon cavity. 

The response of the system described through the linearized optomechanical Hamiltonian could be alternatively 
obtained following the quantum noise approach (see the discussion in the manuscript).^ While this approach is usually 
applied to derive the limit for the phonon population in the cooling setup,!^ it can also be extended to describe the 
effect of the heating of the mechanical degree of freedom. 


1. System 

Consider the Hamiltonian given by Eq. (1) in the main text in the frame rotating with the frequency of the laser 
uji (A = ujc — oji) and with h = 1: 


H = Ad^a -I- uJmVh + iH(a^ — a) — go)a(h + b^), (Bl) 

With the inclusion of noise fields Oin, bin, which are taken to have zero mean value and arise from d-correlated thermal 
baths,!^ the dynamics equations for operators d and b are: 

d = —{k/ 2 + iA)d -f fl -I- igd(b + P) + ^/Kain{t), (B2) 

b = -(7m/2 + iuJm)b -I- igd^d + y/^bin{t). (B3) 


2. Linearization procedure 


The Raman scattering setup discussed in this manuscript is characterized by uncommonly large values of the 
optomechanical coupling parameter. In particular, the granularity parameter defined as the ra tio q/ n is of the order 
of 10“^ - such values have so far only been registered in experiments with ultracoldatoms.^^illQl The theoretical 
analysis of the regime where the coupling parameter g and cavity width k are similail^^J^ shows the breakdown of 
this linearization. To relate this result to our work, in Fig. SI we have analyzed this approximation for the value of 
the coupling coefficient used through most of the manuscript ((a), g = go = k/ 250), and significantly increased in (b) 
{g = lOgo) and {c){g = 30go). The strengths of th e Stokes and anti-Stokes signals were obtained by solving the master 
equation with the full (empty squares, Eq. (Bl)) or linearized (full circles, Eq. (lO)in the main text) Hamiltonian. 


These results indicate that in the regime of parameters discussed in the manuscript both Hamiltonians give identical 
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strengths of the Stokes and anti-Stokes scattering. However, with increasing coupling parameter (by a factor of 10 in 
Fig. Sl(b) or 30 in Fig. Sl(c)), the linearized Hamiltonian begins to over-estimate the inelastic scattering. 


(a) 



[eV^] 


• Stokes - num. linearized 

• anti-Stokes - num. linearized 


(b) 



SJ2 [eV^] 

n Stokes - num. full 
□ anti-Stokes - num. full 


(c) 



— Stokes - analytical 

— anti-Stokes - analytical 


FIG. 5. Dependence of the Stokes (blue) and anti-Stokes (orange) scattering on the excitation power and coupling parameter 
g, calculated numerically using full (empty squares, Eq. (Bll) or linearized (full circles, Eq. (10)) Hamiltonian or analytically 


(solid lines) from the quantum regression theorem (QRT) using Eqs.(B20l and (B21| for Stokes and anti-Stokes, respectively. 

iroughout the manuscript go = k/ 250: (a) g = go, (b) 


Coupling coefficients are given as multiples of the coefficient used t 
g = lOgo and (c) g = 30go- The dashed lines denote the linear dependence of the Stokes signal on H 
classical theories, obtained by taking the first term of Eq. (B201 exclusively. 


expected from the 


In the Methods section, we have introduced the linearized form of the Hamiltonian which does not include any 
coherent driving terms of the cavity or the vibrations. To this goal, we have defined respective coherent amplitudes 
a's and /3' (see Eqs. (12) and (13)). Inserting Eq. (12) into (13) and noting that only the real part Re(/3g) is relevant, 
we arrive at a cubic equation, namely 


V[Re(/3()]3_4gA[Re(/3')]^ 


(f)‘ 




MP's) - 


guj, 




iim/2y 


= 0 . 


(B4) 


The leading order of Eq. (B4) gives the approximations to /?' and a'g, which can be otherwise derived as the coherent 
amplitudes of the undisplaced operators (b) and (a), respectively. The subsequent corrections decrease very quickly 
and are negligible for our parameters. However, the presence of these terms, and in particular the appearance of 
higher powers of is evidence of the fundamental non-linearity of our system. 

Finally, we note that displacing the original Liouvillian by a), instead of (and neglecting the cubic term in the 
Hamiltonian only then) allows our model to also account for higher-order interactions between the plasmons driven 
by the pumping laser and the coherent phonons generated via the cubic interaction. 

In our numerical calculations, we have implemented the Hamiltonians given by Eq. (II) with and without the 
additional non-linear term oc a^a, and solved the corresponding master equation, describing the vibrational and 
plasmonic degrees of freedom in the basis of Fock states. For the strongest pumping and coupling parameters, the 
calculations were ensured to converge by using up to 15 and 10 Fock states for the description of the vibrational and 
photonic state, respectively. 


3. Steady state of the linearized Hamiltonian 


From here on, we will focus our attention on the dynamics of system given by the linearized Hamiltonian H' and 
the master equation (Eq. (2) in the main text). We can thus rewrite the Heisenberg equations given in Eq. (B2) 
using the displaced photonic operator a as 


a = —(k/ 2 -I- fA)d -I- igasib P) ^/nainit), 


(B5) 
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b = + iuJm)b + ig{a*a + a^a''') + ifflasP + V^^in(0- (B6) 

Denoting by A a column vector with expectation values of operators (a^ ,b,,b^, b) and by Fab the correlation matrix 

/at\ 


Fab = 


a 
b^ 

V b j 


^ a b b^ 


(B7) 


we can rewrite Eqs. (B5) and (B6) as 


d 

dt 


where 


A = MA + D, 

(B8) 

'ab = MVab + TabM^ + AD^ + D A^ + E, 

(B9) 

D = i5|asP(0,0,-l,l)^, 

(BIO) 

E = diag([0, K, 7„(1 -f nf)]), 

(Bll) 


for diag denoting a diagonal matrix and M the dynamical matrix 

f-{K/2-iA) 0 -igal -iga*s 

0 —(K/2 + tA) igois igcts 

-igas -igat -( 7^/20 

\ igas iga*s 0 -( 7^/2 + / 


M = 


(B12) 


Denoting by fab = Fab — AA'' the covariance matrix, and putting for simplicity nf^ = 0, we get a simple equation of 
motion; 


—fab = Mfab + fabM'l' + diag([0, K, 0,7^]) ■ 


(B13) 


Computations can be simplified by vectorizing these equations. To this goal, we denote by Fab the column vector 
formed by the rows of fob (and, similarly, by E the column vector formed by the rows of E) and by M the 16 x 16 
matrix M 0 1 + 1 0 M*. Then Eq. (B13) reads 


UTTab — ait ab + E. 
dt 

The displacements and covariances at time t are then given by 

A{t) = -M-^D + e‘^ (A(0) + M-^D) 
fab(t) = -M-^E + (fab(0) + M-^E 

and thus the steady-state displacements and covariances are 

Ass = -M-^D 
Tab,SB = -M-^E. 


(B14) 

(B15) 

(B16) 

(B17) 

(B18) 


4. Quantum regression theorem 

We can obtain correlation functions by noting that, from the QRT, the equation of motion for {X{t)Y{Q)) is the 
same as that for {X{t)) and we find for the matrix 


C(t) = 


f b-Ht) \ 

d{t) 
6t(t) 

V kt) j 


(a(0) at(0) 6(0) 6l'(0)))ss, 
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where the subscript ss denotes the steady-state, C{t) = MC{t) -|- DA^ and thus 

C{t) = -M-^DAl + e*^^Tab,ss, (B19) 

for the steady-state covariance matrix Fa&.ss = Fab.ss — ^ss^L- 

Up to the constant terms we then find that the two-time correlator (a^(t)o(O)) used to calculate the spectrum of 
emission from the cavity is given by the first element of the vector defined as exp(tM) acting on the first column of the 
steady-state covariance matrix Fab.ss- We note that the time-dependent contribution to C{t) depends solely on Fab,ss, 
which is independent of the displacement ag. Therefore, the time-dependent part of (a^(t)a(0))ss is independent of 
whether we define it through the displaced or un-displaced cavity operators. 

For the resonant case (A = 0) and without thermal pumping (T = 0 K) the exponent of M can be found analytically, 
and we can rewrite the above product as a series of expressions with exponential factors given by the eigenvalues of M: 
{g-K/ 2 ,g-( 7 m/ 2 ±ii,;„)}^ From those, we choose the terms oscillating at frequencies iwm which decay as as 

they govern the strengths of Stokes and anti-Stokes scattering. After some algebra, we obtain the explicit expressions 
for the Stokes and anti-Stokes emission as 3rd order polynomials of 

A . ^ 4 

S{ujs) = ^(52^2 -t S4U^), (B20) 


Ai 

S{uJas) = (B21) 

The exact formulas for Si and are somewhat lengthy, but take a simpler form in the limit 7m ^ k, ujm and are then 
given by 


S2 


S4 


04 



_ 

22COm I 



__ 

2zCUm I 


S4. 


2 


4 

K 9 K 

- = S2-, 

'Jm 7m 


(B22) 

(B23) 

(B24) 


The first term in the expression for S{ujs), proportional to is dominant for low pumping power. For larger U 
the second term becomes dominant, and yields the non-linear dependence discussed in the manuscript. Note that the 
coefficients of the terms proportional to are equal both for Stokes and anti-Stokes, and thus for larger pumping 
powers the Stokes/anti-Stokes ratio becomes asymptotically independent of U. 

These two functions are plotted in Fig. SI with blue (Stokes, S{ujs)) and orange (anti-Stokes, S{ujas)) lines and 
match perfectly with the filled circles representing the numerical solution of the linearized Hamiltonian. 

If T > 0, additional terms appear in the expressions for the Stokes and anti-Stokes emission, reflecting (i) the 
thermal pumping of the molecule, which effectively changes the strength of the anti-Stokes scattering and (ii) the 
enhancement of the Stokes scattering by stimulated phonon emission. Specifically, we find that Eq. (B20) can be 
rewritten 


A, ,4 

S(u;s) = -^ls 2 (l+nl^)n^ + 54^4], (B25) 

and the anti-Stokes emission becomes 

A, .4 

S(u^as) = — —(a2nl^n^ + a4n^), (B26) 

Im 


where 02 = S 2 . Thus, in the region of thermal pumping (where we neglect all the terms oc in the above equations), 
we recover the well-known formula for the anti-Stokes/Stokes raticJ^l 


SjuJas) _ / 

S{ujs) V 


n 


th 

h 


^aS 

UJS 


4 


1 -I- ' 


(B27) 
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5. Phonon population 

Our analytical approach allows us to write down the explicit expression for the incoherent phonon populations 
^mcoh thermal nl^ and vibrational pumping, and coherent phonon population n™** in the case of the laser 

tuned to the cavity resonance (A = 0): 


rib 




^incoh 


+ n' 


coh 


where 


^incoh 


= n\!' + 


/ 4gfl y K + 'lm 
V « / 7m|« + 7m + 


^th 


S2^ 


Im 


(B28) 


(B29) 


o^coh 

rib 



(B30) 


6. Quartic dependence of Stokes intensity on the enhancement of the incident field 


By inserting the definitions of (Eq. (A7)) and g (Eq. (Al)) into the lowest-order expression for S{ujs) (Eq. (B20)), 
we can rewrite S{uis) oc K'^ jVeS, where K is the enhancement of the incid ent field at the position of the molecule. 
The inve rse volume factor can be, by relating the Purcell factor (Eq. (A2)) to through the reciprocity theorem 
(see Eq. @),™ shown to be proportional to K^, allowing us to recover the expected dependence S{ujs) oc AT^. 

It should be also noted that a similar dependence of the anti-Stokes intensity on K can be retrieved only in the 
thermal pumping regime, where the phonons are primarily provided by the thermal bath, and the anti-Stokes intensity 
is proportional to In the vibrational pumping regime, the phonons are provided by the Stokes transitions, and 

thus we expect to retrieve a higher-order dependence of S[uJas) on K. 


7 . Dynamical backaction and cooperativity 

For the laser detuned {A = ujc — uJi 0) from the plasmonic cavity, the incoherent population of phonons - both 
vibrationally and thermally pumped - is further enhanced (for A < 0) or quenched (for A > 0) by the onset of the 
dynamical backaction mechanism.l^ 

In this framework, we have introduced in the text the cooperativity-like parameter C 

(7=^. (B31) 

In the limit of C reaching —1, we observe the onset of phonon lasing^^^ in which the Fock basis required for the 
description of the final phonon state would be infinitely large. However, in the regime of parameters discussed in our 
work, C is far from this value. To illustrate this effect, in Fig. S2 we show the values of cooperativity calculated for 
the pumping intensities used in Fig. 3 in the main text of the manuscript. 


Appendix C: Threshold for the onset of the phonon-stimulated Raman scattering 


From Eq. (3) of the main text of the manuscript, we can derive an approximate criterion for the onset of the 
phonon-stimulated Raman scattering for the laser tuned to the cavity resonance, defined arbitrarily by the stimulated 
Stokes emission becoming as strong as the spontaneous emission: 

nt + S2^^—>l. (Cl) 

Im 


Dropping the first term on the left-hand side, which becomes comparatively very small as the power is increased, and 
inserting the definitions of S 2 , g and H from Eqs. (Al), (B22| and (A9), we can show that 


52^2^— oc 


R.Q[ 

Kff 


riXc 


\k - 2iuirn\^ 


\En 
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FIG. 6. Absolute values of C for the pumping intensities from = 10 ^ eV^ (red lines) to = 0.25 eV^ (green lines). 
All the parameters of the system are identical to those discussed in the main text of the manuscript. 


Therefore, the phonon-stimulated processes can become relevant for weaker intensities of the incident laser power if 
we consider either a molecule with a stronger Raman activity and a weaker decay rate of phonons 7m, or a better 
cavity with smaller effective volume and decay rate k or larger quantum yield rj. 


Appendix D: Two-photon physical spectrum 


An important consequence of considering the full quantum model of the SERS is the possibility of studying quantum 
correlations induced between the phonon and cavity mode. A way of characterizing these correlations may be through 
quantum state tomography. However, experimentally accessing the intracavity or the phonon field is impossible with 
the current technology. 

An interesting alternative consist of studying the correlations between the Stokes/anti-Stokes photons in the outfield 
which can be done by inserting filters into the st anda rd Hanbury-Brown-Twiss setujJ ^. Recen tly, it has been shown 
how the correlations between different frequencie^^^EHEIl or their full landscapes (wi, a; 2 )P^^ 2 ^shed light on quantum 
dynamics that is otherwise hidden in standard spectroscopy or color-blind correlations such as violation of classical 
Bell and Cauchy-Schwartz inequalities.^^ Moreover, recent advances with streak camera setups allows to do such 
correlations in the ps scale.l^ 

Interestii^y, the regions of strongly bunched frequency correlations have been linked to the violation of classical 
inequalities^ and can be optimized through proper filter engineering^, potentially providing access to having a 
quantum correlated emission from our setup. 

Theoretically, these time and frequency-resolved photon correlations are computed through intensity-intensity cor¬ 
relations: 


= lim 




t—>-oo 


((A 




(Dl) 


where is the output field after passing through a Lorentzian frequency filter 

with central frequency oji, and width T^, at time t. We note that throughout this section d denotes the original 
operators of the cavity, before displacement. In principle, in order to compute these correlations one must apply 
the quantum regression theorem three times and perform the integrals afterwards. However, in this work we use the 
method recently developed in Ref. 1501 that avoids this complication by coupling the mode of interest, i.e., a in our 
case, to two-level systems that will play the role of sensors (see the schematic in Fig. S3), with frequencies uJi and 
lifetimes T^, through the following Hamiltonian: 

^sens = ^ + d(T^^), (D2) 

2 ^ 1,2 


where the coupling must be sufficiently weak so that the dynamics of the sensors does not perturb the system., i.e., 
Aef/Ti <C 7s, where 7 s is the smallest transition rate of interest (here 7 s = 7m)- Notice that the following can always 
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be imposed as is a free non-zero parameter that we can choose at will in our simulation. Under this assumption it 
can be shown that: 



(: T[ni(0)w2(T)] :)ss 

(ni) SS (^2) SS 


(D3) 


(i) (i) 

with hi = dy (T_ . This simplifies the calculation at the cost of a small increase in the dimension of the Hilbert 
space. For example, in the case of coincidences, i.e., t = 0, that we are mainly interested in, we only need to compute 
one-time correlators, avoiding the need of the quantum regression theorem. 






FIG. 7. Two-photon frequency-resolved correlators gp^p(cu 5 , tUas) are calculated by 
(eSiCaS <C VTiiiF/ 2) two-level sensors and calculating their intensity correlations.® 


adding two lossy (F) weakly coupled 


Note that the state of the Stokes/anti-Stokes output modes expected in our setting is to good approximation given 
by a very weak two-mode squeezed vacuum state {ip) « (1 -|- edsO-ls) |0) reflecting the fact that every anti-Stokes 
photon must (at T = 0 K) be preceded by a Stokes photon. With the use of alone this quantum-correlated 
(entangled) state is, however, not distinguishable from the classically correlated mixture (1 — e^) |0) (0| -I- |2) (2|, 

both of which show diverging as e —>■ 0+. 


1. Dependence of Stokes-anti-Stokes correlations on the parameters of the system 


As we show in Fig. 4 in the manuscript, the Stokes-anti-Stokes correlations decrease both with increasing tem¬ 
perature and increasing pumping intensity. Furthermore, the calculated correlations are significantly weaker for the 
smaller coupling parameter g. 

To qualitatively describe these dependencies, we consider the Stokes-anti-Stokes correlation (wg, WaS, "t = 0) 
defined by the joint probability of detecting a pair of photons, simultaneously reaching the Stokes (SD, see schematics 
in Fig. S4) and anti-Stokes detectors (aSD), normalized by the probabilities of the SD and aSD detections 


9rj'(^s,^aS,T = 0 ) 


Pis, aS) 
PiS)PiaS)' 


(D4) 


Recently, Kasperczyk et aZ.l^and Parra-Murillo et al.^^ha-ve analyzed an idealized version of such detection setup, in 
which the detection probabilities PiS) and P(aS') are proportional to the respective intensities of the Raman intensities 
Siojs) and Siojas), respectively. They showed that if all anti-Stokes photons are produced as part of correlated pairs 
of Stokes and anti-Stokes photons emitted from the cavity (as represented pictorially in Fig. S4(a)), then the joint 
probability PiS, aS) should be proportional to the intensity of the anti-Stokes emission PiS, aS) oc P(aS'), resulting in 
the correlation g^"^^ inversely proportional to the probability of the Stokes emission (g^^) oc 1/P(S') oc I/Siujs) oc 
where the last proportionality is due to first order processes such as the one depicted in Fig. S4(b)). 
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FIG. 8. Main processes contributing to the correlation 5 r^^^(ais, oJaS, T = 0) between Stokes and anti-Stokes photons emitted from 
a SERS cavity observed by detectors of Stokes (SD) and anti-Stokes photons (aSD). For weak pumping and zero temperature, 
the dominant process creating anti-Stokes photons is (a) the correlated emission of Stokes and anti-Stokes photon generated by 
exchanging a single phonon, while (b) Stokes photons are mainly created by emitting a phonon that decays without producing 
an anti-Stokes transition, (c) At non-zero temperatures, the anti-Stokes photons reaching the detector can also be generated 
through the absorption of thermal phonons, (d-e) The detectors characterized by a non-zero spectral width F can also detect 
elastically scattered photons instead of the Stokes or anti-Stokes photons. 


In a realistic scheme, as discussed in the manuscript, this simple relationship is modified by the presence of additional 
process es (n otably those illustrated in Fig. S4(c)-(e)) that change the dependence on of the probabilities appearing 
in Eq. (D4). 


1. The anti-Stokes detector (aSD) observes photons emitted due to the thermal population of phonons (T > 0), 
(see Fig. S4(c)). This effect is particularly important for low pumping intensity, where the thermal contribution 
dominates the phonon population. 


2. The Lorentzian filters used in our detection scheme do not perfectly select the Stokes and anti-Stokes photons, 
breaking the assumption that P{S) oc S{u}s) and P{aS) oc S{uJas) (see Fig. S4(d) and (e)). The importance 
of this effect can be reduced by using different type of detectors, cut-off filters or narrower Lorentzian filters 
(smaller F). 

For sufficiently weak pumping intensity these additional processes modify the dependence of In particular, 
both processes add a oc 17^ contribution to P{aS) and lead to a plateau of the correlations in the weak-pumping 
regime (see Fig. 4(a) in the main text). 

Finally, we would like to point that the correlations of emitted light could be potentially modified by the presence 
of a second, overlapping cavity mode. A similar effect has been reported recently in the study of a bi-modal cavity 
coupled to a two-level system .1^ 
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